from IPython.display import HTML
HTML('''<script>
code_show=true;
function code_toggle() {
if (code_show){
$('div.input').hide();
} else {
$('div.input').show();
}
code_show = !code_show
}
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Mostra codice/ Nascondi codice."></form>''')
from sklearn.model_selection import train_test_split
from sklearn.neighbors import KNeighborsClassifier
from sklearn import metrics
import pydotplus
import matplotlib.colors as mcolors
import statsmodels.formula.api as smf
from sklearn.metrics import mean_squared_error
from IPython.core.display import HTML
from sklearn.model_selection import train_test_split, cross_val_score, cross_val_predict, StratifiedKFold
import plotly.express as px
from sklearn.model_selection import train_test_split, GridSearchCV, cross_val_score, RepeatedStratifiedKFold, StratifiedKFold
from sklearn.metrics import accuracy_score, confusion_matrix,roc_curve, roc_auc_score, precision_score, recall_score, precision_recall_curve
from sklearn.metrics import f1_score
import warnings
warnings.filterwarnings('ignore')
from sklearn.metrics import fbeta_score
from sklearn.metrics import f1_score
from sklearn.metrics import precision_score
from sklearn.metrics import recall_score
from sklearn.svm import SVC
from sklearn.preprocessing import OneHotEncoder
from sklearn.model_selection import StratifiedShuffleSplit
from sklearn import preprocessing # to normalise existing X
from sklearn.cluster import KMeans
from sklearn.preprocessing import LabelEncoder
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn import preprocessing as prp
import plotly.express as px
import ipywidgets as widgets
import plotly.graph_objs as go
from sklearn.linear_model import LinearRegression
import statsmodels.formula.api as stat
from sklearn.model_selection import train_test_split
lm =LinearRegression()
from plotly.offline import plot, iplot
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import classification_report
from sklearn.metrics import confusion_matrix,accuracy_score
Analisi descrittiva¶
data = pd.read_csv('C:/Users/Simona/OneDrive/Desktop/MAGISTRALE/Secondo semestre/Business Analytc/casi/pilgA/PilgrimABC.csv', sep=',')
Controlliamo le dimensioni del dataset e verifichiamo la presenza di valori mancanti nei dati.
print('osservazioni, variabili =', data.shape[0], 'x', data.shape[1])
print('Dati mancanti per variabile')
print(data.isnull().sum())
Decidiamo di eliminare le colonne '0Billpay', '9Billpay' per la nostra analisi e introduciamo alcune nuove variabili al fine di gestire i valori mancanti per gli attributi '9Age', '9Inc', '0Profit' e '0Online'.
In particolare introduciamo le seguenti nuove variabili:
- 'AgeGiven', che vale 1 se l'età del cliente è nota;
- 'Retain', che vale 1 se il cliente ha abbandonato la banca;
- ' AgeZero', che assume lo stesso valore di '9Age' per i dati presenti e 0 per quelli mancanti;
- 'IncomeGiven', che vale 1 se è noto il reddito del cliente:
- 'Profit0', che assegna 0 ai dati mancanti di 'Profit0';
- 'IncomeZero', che assegna 0 ai dati mancanti in '9Inc';
- 'AgeAverage' , che assume lo stesso valore di '9Age' per i dati presenti e la media degli altri valori per quelli mancanti.
data = data.drop(columns=['0Billpay', '9Billpay'])
data['AgeGiven']= np.where(data['9Age'].isnull(), 0, 1) #1 ho l'informazione della etÃ
data['Retain']= np.where(data['0Profit'].isnull(), 1, 0) #0 ho conservato il cliente
data['AgeZero'] = data['9Age'].fillna(0)
data['IncomeGiven']= np.where(data['9Inc'].isnull(), 0, 1) #1 ho l'informazione della data
data['Profit0'] = data['0Profit'].fillna(0)
data['IncomeZero'] = data['9Inc'].fillna(0)
data['AgeAverage'] = data['9Age'].fillna(data['9Age'].mean())
data['Online0'] = data['0Online'].fillna(0)
data.rename(columns={"9Profit": "Profit9", "9Online": "Online9", "9Tenure": "Tenure9", "9District":"District9"}, inplace=True)
#data.describe()
data2= data.copy()
data2= data2.drop(['ID','9Age','9Inc','0Online', '0Profit'], axis=1)
data=data2.copy()
Ecco il dataset, dopo le modifiche:
data
Codifichiamo le variabili categoriche e verifichiamo che il dataset sia bilanciato, rispetto alla variabile 'Retain'.
columns_to_encode = ['District9','AgeZero','IncomeZero']
encoder=OneHotEncoder(sparse=False)
data_pre= data
X_encoded = pd.DataFrame (encoder.fit_transform(data_pre[['District9']]))
X_encoded.columns = encoder.get_feature_names(['District9'])
data_pre.drop(['District9'] ,axis=1, inplace=True)
p1= pd.concat([data_pre, X_encoded ], axis=1)
X_encoded = pd.DataFrame (encoder.fit_transform(p1[['AgeZero']]))
X_encoded.columns = encoder.get_feature_names(['AgeZero'])
p1.drop(['AgeZero'] ,axis=1, inplace=True)
p2= pd.concat([p1, X_encoded ], axis=1)
X_encoded = pd.DataFrame (encoder.fit_transform(p2[['IncomeZero']]))
X_encoded.columns = encoder.get_feature_names(['IncomeZero'])
p2.drop(['IncomeZero'] ,axis=1, inplace=True)
data_enc= pd.concat([p2, X_encoded ], axis=1)
data_enc.rename(columns={'AgeZero_1.0':'AgeZero_1fascia',
'AgeZero_2.0':'AgeZero_2fascia',
'AgeZero_3.0':'AgeZero_3fascia',
'AgeZero_0.0':'AgeZero_0fascia',
'AgeZero_4.0':'AgeZero_4fascia',
'AgeZero_5.0':'AgeZero_5fascia',
'AgeZero_6.0':'AgeZero_6fascia',
'AgeZero_7.0':'AgeZero_7fascia'},
inplace=True)
data_enc.rename(columns={'IncomeZero_1.0':'IncomeZero_1fascia',
'IncomeZero_2.0':'IncomeZero_2fascia',
'IncomeZero_3.0':'IncomeZero_3fascia',
'IncomeZero_0.0':'IncomeZero_0fascia',
'IncomeZero_4.0':'IncomeZero_4fascia',
'IncomeZero_5.0':'IncomeZero_5fascia',
'IncomeZero_6.0':'IncomeZero_6fascia',
'IncomeZero_7.0':'IncomeZero_7fascia',
'IncomeZero_8.0':'IncomeZero_8fascia',
'IncomeZero_9.0':'IncomeZero_9fascia'},
inplace=True)
data_enc.rename(columns={'District9_1300.0':'District_thir',
'District9_1200.0':'District_sec',
'District9_1100.0':'District_first'},
inplace=True)
data_clas= data_enc.drop(columns=['Online0'])
group_size=data['Retain'].value_counts().sort_index()
group_names=['Clienti Fedeli (' + '{:1.1f}'.format(group_size[0]/(group_size[0]+group_size[1])*100) + '%)',
'Clienti Persi (' + '{:1.1f}'.format(group_size[1]/(group_size[0]+group_size[1])*100) + '%)']
a, b=[plt.cm.Greens, plt.cm.Reds]
fig, ax = plt.subplots()
ax.axis('equal')
mypie, _ = ax.pie(group_size, labels=group_names, colors=[a(0.6), b(0.6)], explode=(0.1,0),
textprops={'size': 18, 'color':'black'})
plt.figsize=(5,5)
plt.setp( mypie, width=1, edgecolor='white')
plt.margins(0,0)
fig = plt.gcf()
fig.set_size_inches(4,4)
plt.show()
Come si nota il dataset è sbilanciato, questo implica la necessità di utilizzare, in fase di classificazione, determinate tecniche che tengano conto della distribuzione del dataset.
Visualizziamo il dataset mediante alcuni grafici, ad esempio la distribuzione delle diverse fasce d'età nel presenti nel dataset.
a=31634
#vediamo se è bilanciato
group_size=data['AgeAverage'].value_counts().sort_index()
group_names=['Fascia-età 1 (' + '{:1.1f}'.format(group_size[1]/(a)*100) + '%)',
'Fascia-età 2 (' + '{:1.1f}'.format(group_size[2]/(a)*100) + '%)',
'Fascia-età 3 (' + '{:1.1f}'.format(group_size[3]/(a)*100) + '%)',
'Fascia-età 4 (' + '{:1.1f}'.format(group_size[4]/(a)*100) + '%)',
'Età assente (' + '{:1.1f}'.format(8289/(a)*100) + '%)',
'Fascia-età 5 (' + '{:1.1f}'.format(group_size[5]/(a)*100) + '%)',
'Fascia-età 6 (' + '{:1.1f}'.format(group_size[6]/(a)*100) + '%)',
'Fascia-età 7 (' + '{:1.1f}'.format(group_size[7]/(a)*100) + '%)']
#colorlist = ['crimson','k','r','gold','g','b','c','d','a','l']
colors = mcolors.TABLEAU_COLORS.values()
fig, ax = plt.subplots()
ax.axis('equal')
mypie, _ = ax.pie( group_size, radius=2.5, labels=group_names, colors=colors,
textprops={'size': 18, 'color':'black'})
plt.setp( mypie, width=0.7, edgecolor='white')
plt.margins(0,0)
fig.set_size_inches(3,3)
plt.show()
#grafic
#grafico2
x1 = data_enc.loc[data_enc.AgeGiven==0, 'Profit9']
x2 = data_enc.loc[data_enc.AgeGiven==1, 'Profit9']
fig, axes = plt.subplots(1, 2, figsize=(10, 3), sharey=True, dpi=100)
axes[0].title.set_text('Dato su età mancante')
axes[1].title.set_text('Dato su età presente')
sns.distplot(x1 , color="dodgerblue", ax=axes[0], axlabel='Profitto9')
sns.distplot(x2 , color="deeppink", ax=axes[1], axlabel='Profitto9')
#plt.hist(x1, alpha=0.5, bins=100, density=True, stacked=True, label='age', color='dodgerblue')
#plt.hist(x2, alpha=0.5, bins=100, density=True,stacked=True, label='noage', color='deeppink')
plt.tight_layout();
barWidth = 0.25
plt.figure(figsize=(10,5))
# set width of bar
barWidth = 0.4
a= data_enc.loc[data_enc.IncomeGiven==0, 'AgeAverage']
b= data_enc.loc[data_enc.IncomeGiven==1, 'AgeAverage']
a= group_size=a.value_counts().sort_index()
b= group_size=b.value_counts().sort_index()
# set height of bar
bars1 = a
bars2 = b
# Set position of bar on X axis
r1 = np.arange(len(bars1))
r2 = [x + barWidth for x in r1]
# Make the plot
plt.bar(r1, bars1, color='g', width=barWidth, edgecolor='white', label=' INCOME mancante')
plt.bar(r2, bars2, color='orange', width=barWidth, edgecolor='white', label=' INCOME presente')
# Add xticks on the middle of the group bars
plt.xticks([r + barWidth for r in range(len(bars1))], ['Age1', 'Age2', 'Age3', 'Age4', 'Età mancante','Age5','Age6','Age7'], fontweight='bold', fontsize=20)
plt.xticks(rotation=90)
plt.yticks(fontsize=20)
# Create legend & Show graphic
plt.legend(fontsize=15)
plt.show()
Il grafico precedente mostra che nella maggior parte delle osservazioni in cui non è presente il dato su Income, non è presente neppure il dato su Età .
sns.set(rc={'figure.figsize':(15,10)})
sns.set(font_scale = 2)
bplot=sns.boxplot(y='Profit9', x='Retain',
data=data_enc,
width=0.5,
palette="nipy_spectral")
Il boxplot mostra, in entrambe le classi ( Retain=1 o Retain=0) , la presenza di numerosi (eventuali) outliers .
data= data2.copy()
Regressione¶
Cerchiamo di capire, grazie ad alcune semplici regressioni lineari, in che modo le variabili 'Age' , 'Income' e 'Online'influenzino il profitto .
Ruolo della variabile 'Age'¶
Decidiamo come trattare la variabile Age. Abbiamo diverse alternative:
- come variabile categorica, a cui associamo la fascia 0 ai dati in cui l'età è mancante
- come variabile numerica:
- in cui la media è sostituita ai dati mancanti ;
- in cui 0 è sostituito ai dati mancanti .
X= data.drop(['Retain', 'Profit9','Profit0'], axis=1)
Y= data['Profit9']
data_reg= X.join(Y)
Proviamo a vedere la relazione tra 'Profit9' e 'AgeGiven'.
modello1 = smf.ols('Profit9 ~ AgeGiven', data=data_enc).fit()
def short_summary(est):
return HTML(est.summary().tables[1].as_html())
print('Il valore di MSE del modello vale' + ' {:1.2f}'.format(modello1.mse_resid)+ ', mentre l\'R^2adj vale' + ' {:1.4f}'.format(modello1.rsquared_adj))
short_summary(modello1)
Il modello precedente mostra che la variabile 'AgeGiven è significativa (p-value < 0.05) . In particolare il coefficiente è positivo per cui ci aspettiamo che il profitto per i dati di cui è nota l'età sia superiore, in media , rispetto al profitto delle osservazioni di cui non è nota l'età .
Infatti se vogliamo entrare più in dettaglio, il valore di intercetta non è altro che la media della variabile 'Profit9' per i dati in cui 'AgeGiven' è uguale a 0.
a= data_enc[data_enc['AgeGiven']==0]
a["Profit9"].mean()
Invece, il coefficiente $\beta$ di AgeGiven non è altro che la differenza tra l'intercetta e la media della variabile 'Profit9' per i dati in cui 'AgeGiven' è uguale a 1.
b= data_enc[data_enc['AgeGiven']==1]
b["Profit9"].mean()- a["Profit9"].mean()
Adesso, proviamo a vedere la relazione tra 'Profit9' e la variabile 'Age' come categorica.
modello1 = smf.ols('Profit9 ~ AgeZero_1fascia+ AgeZero_2fascia+ AgeZero_3fascia+ AgeZero_4fascia+ AgeZero_5fascia+ AgeZero_6fascia + AgeZero_7fascia', data=data_enc).fit()
print('Il valore di MSE del modello vale' + ' {:1.2f}'.format(modello1.mse_resid)+ ', mentre l\'R^2adj vale' + ' {:1.4f}'.format(modello1.rsquared_adj))
short_summary(modello1)
#rint(results.summary())
#Y_pred= modello1.predict(X)
#mse = mean_squared_error(Y, Y_pred)
#_pred= results.predict(X_test)
#mse = mean_squared_error(Y_test, Y_pred)
#results = smf.ols('Profit9 ~ Online9+ AgeGiven', data=data_reg).fit()
#print(results.summary())
Come si nota il valore dell'intercetta di questo modello è uguale al valore dell'intercetta del modello precedente. Questo non ci stupisce, perchè in questo caso la categoria di riferimento della regressione è proprio la classe di clienti di cui non è nota l'età . Dai coefficienti della regressione abbiamo una informazione in più sulle medie della variabile 'Profit9'. Infatti i coefficienti sono negativi per le categorie 1 e 2, il chè implica che la media delle variabile 'Profit9' è minore rispetto alla categoria di riferimento.
Mentre i clienti delle restanti categorie sono, in media, più profittevoli rispetto a coloro di cui non è nota l'età . A titolo esemplificativo, vediamo il calcolo del coefficiente relativo a AgeZero_1fascia:
b= data_enc[data_enc['AgeZero_1fascia']==1]
b["Profit9"].mean()- a["Profit9"].mean()
Adesso proviamo a considerare la variabile 'AgeAverage' ( numerica) .
modello2 = smf.ols('Profit9 ~ AgeAverage', data=data_reg).fit()
print('Il valore di MSE del modello vale' + ' {:1.2f}'.format(modello2.mse_resid)+ ', mentre l\'R^2adj vale' + ' {:1.4f}'.format(modello2.rsquared_adj))
short_summary(modello2)
Infine consideriamo la variabile 'AgeZero' (numerica) .
modello3 = smf.ols('Profit9 ~ AgeZero', data=data_reg).fit()
print('Il valore di MSE del modello vale' + ' {:1.2f}'.format(modello3.mse_resid)+ ', mentre l\'R^2adj vale' + ' {:1.4f}'.format(modello3.rsquared_adj))
short_summary(modello3)
Scegliamo di trattare la variabile sull'età come categorica sia perchè sembra essere la scelta più ragionevole, confermata anche dal valore di MSE e R^2adj dei modelli precedenti.
Ruolo della variabile 'Income'¶
Vediamo la relazione tra 'Profit9' e 'IncomeGiven'.
modello1 = smf.ols('Profit9 ~ IncomeGiven', data=data_enc).fit()
print('Il valore di MSE del modello vale' + ' {:1.2f}'.format(modello1.mse_resid)+ ', mentre l\'R^2adj vale' + ' {:1.4f}'.format(modello1.rsquared_adj))
short_summary(modello1)
La variabile 'IncomeGiven' sembra essere significativa ed, in particolare, il profitto è, in media , maggiore per le osservazioni in cui è presente il dato su 'Income'.
Adesso, proviamo a vedere la relazione tra 'Profit9' e la variabile 'Income' come categorica.
modello1 = smf.ols('Profit9 ~ IncomeZero_1fascia+ IncomeZero_2fascia+ IncomeZero_3fascia+ IncomeZero_4fascia + IncomeZero_5fascia+ IncomeZero_6fascia+ IncomeZero_7fascia+ IncomeZero_8fascia+ IncomeZero_9fascia', data=data_enc).fit()
print('Il valore di MSE del modello vale' + ' {:1.2f}'.format(modello1.mse_resid)+ ', mentre l\'R^2adj vale' + ' {:1.4f}'.format(modello1.rsquared_adj))
short_summary(modello1)
Guardando i coefficienti notiamo che essi sono crescenti, fatta eccezione per la terza fascia. Si potrebbe pensare che il valore di Profitto è, in qualche modo, proporzionale al valore di Income.
Ruolo della variabile 'Online'¶
Cerchiamo di capire, in che modo, il fatto che il cliente utilizzi o meno il servizio Online, incida sul profitto.
#Grafico
on= data[data['Online9']==1]
off= data[data['Online9']==0]
plt.figure(figsize=(15,5))
ax = sns.regplot(x="Profit9", y="Profit0", data=on)
plt.title('Utenti che usano i servizi Online', fontsize=25)
plt.xlabel('Profit_1999', fontsize=20)
plt.ylabel('Profit_2000', fontsize=2)
plt.tick_params( labelsize=15)
plt.figure(figsize=(15,5))
#plt.plot(off['Profit9'], off['Profit0'], marker="o", color='green', linestyle="")
ax = sns.regplot(x="Profit9", y="Profit0", data=off, color='green')
#Grafico
plt.title('Utenti che non usano i servizi Online', fontsize=25)
plt.xlabel('Profit_1999', fontsize=20)
plt.ylabel('Profit_2000', fontsize=20)
plt.tick_params( labelsize=15)
Vediamo la relazione tra 'Profit9' e 'Online9'.
modello2 = smf.ols('Profit9 ~ Online9 ', data=data_reg).fit()
print('Il valore di MSE del modello vale' + ' {:1.2f}'.format(modello2.mse_resid)+ ', mentre l\'R^2adj vale' + ' {:1.4f}'.format(modello2.rsquared_adj))
short_summary(modello2)
Vediamo la relazione tra 'Profit0' e 'Online0' ( per i clienti che non hanno abbandonato la banca).
modello2 = smf.ols('Profit0 ~ Online0', data=data).fit()
print('Il valore di MSE del modello vale' + ' {:1.2f}'.format(modello2.mse_resid)+ ', mentre l\'R^2adj vale' + ' {:1.4f}'.format(modello2.rsquared_adj))
short_summary(modello2)
Per entrambi gli anni, il profitto è, in media, maggiore per i clienti che usano il servizio online.
Inoltre cerchiamo di capire come cambia, in media, il profitto, se un cliente passa da Online9=0 a Online0=1, cioè per i clienti che iniziano a usufrire del servizio online. Per tale analisi, considereremo, per ovvie ragione, solo i clienti che non hanno abbandonato la banca.
L'obiettivo è capire se la variazione di profitto di un cliente sia attribuibile, in parte, al fatto che il cliente abbia iniziato ad usare il servizio Online. Inoltre introduciamo una nuova colonna, che non è altro che la differenza tra i valori di Profit0 e Profit9.
Definiamo anche una nuova variabile ('start') come $ Online0 \cdot (1- Online9) $. Il valore di questa variabile assume 1 solo per i clienti che iniziano ad usare il servizio Online al secondo anno. Infatti :
- Online9=0 , Online0=0 --> start= 0
- Online9=0 , Online0=1 --> start =1
- Online9=1 , Online0=0--> start =0
- Online9=1 , Online0=1 --> start =0
Riassumendo i clienti con 'start'=1 sono i clienti che nel 1999 non usano il servizio Online e che nel 2000 lo usano.
data_online= data[data['Retain']==0]
data_online['diff'] = data_online['Profit0'] - data_online['Profit9']
data_online['start'] = data_online['Online0'] *(1- data_online['Online9'])
data_online
Per i dati che soddisfano ' start= 1' la media della differenza tra i profitti vale:
a= data_online[data_online['start']==1]
a['diff'].mean()
plt.figure(figsize=(13,7))
plt.plot(a['Profit9'], a['Profit0'], marker="o", color='brown', linestyle="")
plt.title('Utenti che hanno iniziato ad usare i servizi Online', fontsize=25)
plt.xlabel('Profit_1999', fontsize=20)
plt.ylabel('Profit_2000', fontsize=20)
plt.tick_params( labelsize=15)
Per i dati che soddisfano ' start= 0' la media della differenza tra i profitti vale:
a= data_online[data_online['start']==0]
a['diff'].mean()
plt.figure(figsize=(13,7))
plt.plot(a['Profit9'], a['Profit0'], marker="o", color='green', linestyle="")
plt.title('Utenti con start=0', fontsize=25)
plt.xlabel('Profit_1999', fontsize=20)
plt.ylabel('Profit_2000', fontsize=20)
plt.tick_params( labelsize=15)
Vogliamo stabilire se la differenza delle due medie è staticamente significativa. Per farlo, utlizziamo il seguente modello: $$ diff= Profit_0- Profit_9 = \beta_0 + \beta \cdot start $$
modello2 = smf.ols('diff ~ start', data=data_online).fit()
print('Il valore di MSE del modello vale' + ' {:1.2f}'.format(modello2.mse_resid)+ ', mentre l\'R^2adj vale' + ' {:1.4f}'.format(modello2.rsquared_adj))
short_summary(modello2)
Vediamo che l'intercetta è la media della differenza tra i profitti dei due anni per i clienti con 'start=0' ( che abbiamo trovato precedentemente). Poichè il p-value del coefficiente 'start' vale 0.002, possiamo affermare, che la differenza tra le due medie è significativa e cioè rifiutiamo l'ipotesi che il coefficiente di 'start' sia zero.
Inoltre, dato che il coefficiente di 'start' è positivo , si potrebbe pensare che risulta conveniente per la banca che il cliente inizi ad usare il servizio online.
Verifichiamo cosa succede per i clienti che hanno smesso di usare il servizio online ( e che non hanno abbandonato la banca). Introduciamo la variabile 'stop'= $ Online9 \cdot (1- Online0) $. Il valore di questa variabile assume 1 solo per i clienti che hanno smesso di usare il servizio Online nel 2000. Infatti :
- Online9=0 , Online0=0 --> stop= 0
- Online9=0 , Online0=1 --> stop = 0
- Online9=1 , Online0=0--> stop = 1
- Online9=1 , Online0=1 --> stop =0
Quello che mi aspetterei è vedere il coefficiente della regressione negativo.
data_online['stop'] = data_online['Online9'] *(1- data_online['Online0'])
Calcoliamo la media della differenza tra i profitti per i clienti che hanno smesso di usare il servizio Online.
a= data_online[data_online['stop']==1]
a['diff'].mean()
Calcoliamo la media della differenza tra i profitti per i clienti con 'Stop=0'.
a= data_online[data_online['stop']==0]
a['diff'].mean()
Contro le mie aspettative, la media della differenza tra i profitti è maggiore per coloro che hanno abbandonato il servizio Online. Ma vediamo se la differenza tra le medie è statisticamente significativa.
modello2 = smf.ols('diff ~ stop', data=data_online).fit()
print('Il valore di MSE del modello vale' + ' {:1.2f}'.format(modello2.mse_resid)+ ', mentre l\'R^2adj vale' + ' {:1.4f}'.format(modello2.rsquared_adj))
short_summary(modello2)
Quello che vediamo è che il coefficiente della regressione è positivo, questo indica che la media per i clienti che hanno smesso di usare il servizio online è più alta ( come abbiamo verificato sopra).
Tuttavia, in questo caso, il p-value del coefficiente 'stop' vale 0.207 > 0.05 e non possiamo rifiutare l'ipotesi che il coefficiente sia zero..
b= data_online[data_online['stop']==1]
plt.figure(figsize=(13,7))
plt.plot(b['Profit9'], b['Profit0'], marker="o", color='brown', linestyle="")
plt.title('Utenti che hanno smesso di usare i servizi Online', fontsize=25)
plt.xlabel('Profit_1999', fontsize=20)
plt.ylabel('Profit_2000', fontsize=20)
plt.tick_params( labelsize=15)
Avendo a disposizione un grande numero di dati, per trovare gli eventuali outliers, calcoliamo gli Z-score delle variabili numeriche e rimuoviamo le osservazioni che hanno un valore di Z-score, in valore assoluto, maggiore di 4.
In particolare, dato che per la regressione di 'Profit9' non utilizziamo come regressore la variabile 'Profit0' ( lo scopo è predire il valore di Profit9 quindi assumiamo di non conoscere il valore di profitto di un anno successivo), le uniche variabili numeriche di cui calcoliamo lo Z-score sono 'Profit9' e 'Tenure9'. Escludiamo dai regressori per la previsione di Profit9 anche la variabile 'Online0'.
Vediamo, mediante i seguenti istrogrammi, quante sono le variabili con uno Z-score maggiore, in valore assoluto, di 4.
import plotly.graph_objects as go
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
#ciao=data_enc[data_enc['Retain']==0]
ciao= data_enc.copy()
l= ciao.shape
salva= ciao.copy()
scaled = scaler.fit_transform(ciao[['Profit9','Tenure9', 'Profit0']])
scaled = pd.DataFrame(scaled, columns=['Profit9','Tenure9', 'Profit0'])
plt.title('Profit9', fontsize=25)
plt.hist(scaled['Profit9'], bins = 50)
plt.xlabel('Zscore', fontweight='bold', fontsize=20)
plt.xticks(np.arange(-2, max(scaled['Profit9'])+1, 1.0), fontsize=20)
plt.yticks(fontsize=20)
plt.figure(figsize=(20,10))
plt.title('Tenure9', fontsize=25)
plt.hist(scaled['Tenure9'], bins = 50)
plt.xlabel('Zscore', fontweight='bold', fontsize=20)
plt.xticks(np.arange(-2, max(scaled['Tenure9'])+1, 1.0), fontsize=20 )
plt.yticks(fontsize=20)
plt.figure(figsize=(20,10))
tol=4
ciao= data_enc.copy()
d = np.logical_or(abs(scaled['Tenure9'])>tol,abs(scaled['Profit9'])>tol)
ciao=pd.concat((ciao,d.rename('col')), axis=1, join='inner')
ciao=ciao[ciao['col']==False]
data_enc=ciao.drop(columns=['col'])
c=data_enc.shape
print("Il numero di osservazioni eliminate è", (l[0]- c[0]))
Dividiamo il dataset in dataset train e test (70%- 30%).
Forward selection¶
Proviamo a predire la variabile 'Profit9'.
data_enc
X= data_enc.drop(['Profit9','Retain', 'Profit0','AgeAverage','Online0'], axis=1)
Y= data_enc['Profit9']
#Y= np.power(Y,2)
X, X_test, Y, Y_test = train_test_split(X, Y, test_size=0.3, random_state=42)
data_reg= X.join(Y)
ciao=[]
n = X.columns.shape[0]
X_train, X_val, Y_train, Y_val = train_test_split(X, Y, test_size=0.5, random_state=42)
selected_features = []
left_feats = X.columns.tolist()
for k in range(1,n+1):
max_r_sq = 0.0
chosen_feat = ""
for feat in left_feats:
model = LinearRegression()
model.fit(X_train[selected_features+[feat]],Y_train)
r_sq = model.score(X_train[selected_features+[feat]],Y_train)
if r_sq > max_r_sq: #if a better model is found (according to R2)
chosen_feat = feat
max_r_sq = r_sq
selected_features.append(chosen_feat)
left_feats.remove(chosen_feat)
print('Le variabili sono selezionate ( in base al valore di R^2) nel seguente ordine:\n')
j=1
for feat in selected_features:
print( '-' + feat)
j +=1
print('\n\n')
best_mse = -50000000000000
best_k = 0
for k in range(1,n+1):
mse = cross_val_score(LinearRegression(), X_val[selected_features[:k]],
Y_val, cv=4, scoring='neg_mean_squared_error')
mse = sum(mse)/5.0 #finding average model mse
if mse>best_mse:
best_mse = mse
best_k = k
print('\n\nIl miglior modello, scelto sulla base del valore di MSE della cross validation (K=5) , è quello con ' + str(best_k) + ' variabili e i coefficienti stimati sono :\n')
model = LinearRegression().fit(X[selected_features[:best_k]],Y)
coefs = model.coef_
Y_pred = model.predict(X[selected_features[:best_k]])
res = Y-Y_pred
q= Y_pred
Y_preds = model.predict(X_test[selected_features[:best_k]])
mse = mean_squared_error(Y_test, Y_preds)
resA= res
yA= Y
for j in range(best_k):
print(selected_features[j].ljust(20,' ') + '{:10.3f}'.format(coefs[j]))
print('\nIl TEST MSE del modello vale ' + '{:1.3f}'.format(mse) + '.')
Come dataset per le seguenti regressioni, consideriamo solo le osservazioni in cui la variabile Retain vale 0, cioè se il cliente non ha abbandonato la banca. Avendo già elimininato gli outliers rispetto alla variabile 'Profit9' e 'Tenure9', non ci resta che eliminarle rispetto a 'Profit0'.
data_enc= data_enc[data_enc['Retain']==0]
data_reg= X.join(Y)
data_enc
l=data_enc.shape
scaler = StandardScaler()
ciao= data_enc.copy()
salva= ciao.copy()
scaled = scaler.fit_transform(ciao[['Profit0']])
scaled = pd.DataFrame(scaled, columns=['Profit0'])
colors = px.colors.qualitative.Plotly
plt.title('Profit0', fontsize=25)
plt.hist(scaled['Profit0'], bins =50)
plt.xlabel('Zscore', fontweight='bold', fontsize=20)
plt.xticks(np.arange(-4, max(scaled['Profit0'])+1, 10) )
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.figure(figsize=(10,7))
tol=4
ciao= data_enc.copy()
a= abs(scaled['Profit0'])>tol
ciao=pd.concat((ciao,a.rename('col')), axis=1, join='inner')
ciao=ciao[ciao['col']==False]
data_enc=ciao.drop(columns=['col'])
c=data_enc.shape
#print("Il numero di osservazioni eliminate è", (l[0]- c[0]))
#c[0]
Dopo aver eliminato gli outliers, nuovamente, dividiamo il dataset in train e test.
X= data_enc.drop(['Profit9','Retain', 'Profit0','AgeAverage'], axis=1)
Y= data_enc['Profit0']
X, X_test, Y, Y_test = train_test_split(X, Y, test_size=0.3, random_state=42)
data_reg= X.join(Y)
Forward selection¶
Esclusione di 'Profit9' dai regressori¶
Proviamo a prevedere Profit0 senza usare Profit9.
ciao=[]
n = X.columns.shape[0] #max number of variables in the model
#cv
#sampling training set and validation set
X_train, X_val, Y_train, Y_val = train_test_split(X, Y, test_size=0.5, random_state=42)
selected_features = [] #empty list where to put selected features
left_feats = X.columns.tolist() #still to be chosen features
for k in range(1,n+1):
max_r_sq = 0.0
chosen_feat = ""
for feat in left_feats: #looking for the next feature to select
model = LinearRegression()
model.fit(X_train[selected_features+[feat]],Y_train)
r_sq = model.score(X_train[selected_features+[feat]],Y_train)
if r_sq > max_r_sq: #if a better model is found (according to R2)
chosen_feat = feat
max_r_sq = r_sq
selected_features.append(chosen_feat) #adding best feature to selected features list
left_feats.remove(chosen_feat) #and removing it from those still to be chosen
print('Le variabili sono selezionate ( in base al valore di R^2) nel seguente ordine:\n')
j=1
for feat in selected_features:
print( '-' + feat)
j +=1
print('\n\n')
#cross-validating and finding the best model
best_mse = -50000000000000
best_k = 0
for k in range(1,n+1):
mse = cross_val_score(LinearRegression(), X_val[selected_features[:k]],
Y_val, cv=4, scoring='neg_mean_squared_error')
mse = sum(mse)/5.0 #finding average model mse
if mse>best_mse:
best_mse = mse
best_k = k
print('\n\nIl miglior modello, scelto sulla base del valore di MSE della cross validation ( K=5) , è quello con ' + str(best_k) + ' variabili e i coefficienti stimati sono:\n')
model = LinearRegression().fit(X[selected_features[:best_k]],Y)
coefs = model.coef_
Y_pred = model.predict(X[selected_features[:best_k]])
#mse = mean_squared_error(Y, Y_pred)
res = Y-Y_pred
Y_preds = model.predict(X_test[selected_features[:best_k]])
mse = mean_squared_error(Y_test, Y_preds)
for j in range(best_k):
print(selected_features[j].ljust(20,' ') + '{:10.3f}'.format(coefs[j]))
print('\nIl TEST MSE del modello vale ' + '{:1.3f}'.format(mse) + '.')
Inclusione di 'Profit9' nei regressori¶
Proviamo a predire profit0 usando profit9
X= data_enc.drop(['Retain', 'Profit0', 'Online0' ,'AgeAverage'], axis=1)
Y= data_enc['Profit0']
X, X_test, Y, Y_test = train_test_split(X, Y, test_size=0.3, random_state=42)
data_reg= X.join(Y)
ciao=[]
n = X.columns.shape[0] #max number of variables in the model
#cv
#sampling training set and validation set
X_train, X_val, Y_train, Y_val = train_test_split(X, Y, test_size=0.5, random_state=42)
selected_features = [] #empty list where to put selected features
left_feats = X.columns.tolist() #still to be chosen features
for k in range(1,n+1):
max_r_sq = 0.0
chosen_feat = ""
for feat in left_feats: #looking for the next feature to select
model = LinearRegression()
model.fit(X_train[selected_features+[feat]],Y_train)
r_sq = model.score(X_train[selected_features+[feat]],Y_train)
if r_sq > max_r_sq: #if a better model is found (according to R2)
chosen_feat = feat
max_r_sq = r_sq
selected_features.append(chosen_feat) #adding best feature to selected features list
left_feats.remove(chosen_feat) #and removing it from those still to be chosen
print('Le variabili sono selezionate ( in base al valore di R^2) nel seguente ordine:\n')
j=1
for feat in selected_features:
print( '-' + feat)
j +=1
print('\n\n')
#cross-validating and finding the best model
best_mse = -50000000000000
best_k = 0
for k in range(1,n+1):
mse = cross_val_score(LinearRegression(), X_val[selected_features[:k]],
Y_val, cv=4, scoring='neg_mean_squared_error')
mse = sum(mse)/5.0 #finding average model mse
if mse>best_mse:
best_mse = mse
best_k = k
print('\n\nIl miglior modello, scelto sulla base del valore di MSE della cross validation ( K=5) , è quello con ' + str(best_k) + ' variabili e i coefficienti stimati sono:\n')
model = LinearRegression().fit(X[selected_features[:best_k]],Y)
coefs = model.coef_
Y_pred = model.predict(X[selected_features[:best_k]])
msetr = mean_squared_error(Y, Y_pred)
res = Y-Y_pred
Y_preds = model.predict(X_test[selected_features[:best_k]])
mse = mean_squared_error(Y_test, Y_preds)
for j in range(best_k):
print(selected_features[j].ljust(20,' ') + '{:10.3f}'.format(coefs[j]))
print('\nIl TEST MSE del modello è ' + '{:1.3f}'.format(mse) + '.')
Osservazioni¶
La presenza di Profit9 come regressore per la previsione di 'Profit0' fa quasi dimezzare il TEST MSE , determinando anche una notevole diminuizione del numero di regressori ( da 11 a 4). Questo fatto non stupisce perchè è chiaro che conoscore il valore del Profitto di un anno è determinante nella previsione del profitto dell'anno successivo. Notiamo, come ci si può aspettare, che la prima variabile selezionata dalla feature selection sia proprio Profit9.
La maggior parte delle variabili selezionate per la previsione 'Profit9' viene usata anche per la previsione 'Profit0' ( primo caso) . Anche questo fatto non stupisce in quanto ci si può aspettare che i fattori determinanti nella previsione del profitto non cambino al trascorrere degli anni.
Analisi dei residui¶
Riportiamo l'analisi dei residui per il primo dei modelli sopra analizzati.
plt.figure(figsize=(10,6))
ax= sns.scatterplot(x=yA, y=resA)
ax.set_xlabel('Profit9 ')
ax.set_ylabel('Residui')
plt.show()
plt.figure(figsize=(10,8))
p= sns.scatterplot(x=q, y=resA)
p.set_xlabel('Valori predetti ( fitted values )')
p.set_ylabel('Residui')
plt.show()
I grafici mostrano che le ipotesi alla base della regressione lineare non sono rispettate. Infatti è evidente un pattern nel grafico dei residui. Ho provato a effettuare alcune trasformazioni delle variabili. In particolare, ho provato ad elevare a diverse potenze le variabili numeriche del dataset, mentre ho escluso la trasformazione logaritmica perchè le variabili assumono anche valore negativo. Queste trasformazioni non hanno tuttavia determinato un miglioramento dei residui.
Di seguito, vediamo i grafici dei residui per la previsione del quadrato di Profit9.
plt.figure(figsize=(10,6))
ax= sns.scatterplot(x=yA, y=resA)
ax.set_xlabel('(Profit9)^2')
ax.set_ylabel('Residui')
plt.show()
plt.figure(figsize=(10,8))
p= sns.scatterplot(x=q, y=resA)
p.set_xlabel('Valori predetti ( fitted values )')
p.set_ylabel('Residui')
plt.show()
Classificazione¶
L'obiettivo è quello di assegnare ai dati una specifica classe ( positiva o negativa) , in base al valore di Retain. Ricordiamo che Retain vale:
- 1 ,se il cliente ha abbandonato la banca -> etichetta positivo ( evento anomalo);
- 0 , se il cliente resta in banca -> etichetta negativo.
In particolare, per la nostra analisi, assumiamo che i due possibili errori di classificazione abbiano un peso diverso. Consideriamo, infatti, che l'errore 'falso negativo' (il cliente abbandona la banca , ma viene classificato come un cliente che resta in banca) abbia una maggiore importanza rispetto a 'falso positivo'. Per sintetizzare, l'obiettivo dei modelli che seguono sarà quello di minimizzare i falsi negativi.
Inoltre, per prima cosa dividiamo il dataset in 'train data' e 'test data'. Quest'ultima porzione di dati verrà utilizzata esclusivamente alla fine della nostra analisi come strumento di confronto tra i diversi algoritmi che analizzeremo. Nello specifico, andremo a dividere il dataset in modo 'stratificato', al fine di avere la stessa distribuzione della variabile ' Retain' nei due dataset che otteniamo.
import itertools
def metrics(true, preds):
"""
Function to calculate evaluation metrics
parameters: true values, predictions
prints accuracy, recall, precision and f1 scores
"""
accuracy = accuracy_score(true, preds)
recall = recall_score(true, preds)
precision = precision_score(true, preds)
f1score = f1_score(true, preds)
f2score= fbeta_score(true, preds, beta=2, pos_label=1)
print ('accuracy: {}, recall: {}, precision: {}, f1-score: {}, f2-score: {}'.format(accuracy, recall, precision, f1score, f2score))
def plot_confusion_matrix(cm,
target_names,
title='Matrice di confusione',
cmap=None,
normalize=True):
accuracy = np.trace(cm) / np.sum(cm).astype('float')
misclass = 1 - accuracy
if cmap is None:
cmap = plt.get_cmap('Blues')
plt.figure(figsize=(8, 6))
plt.imshow(cm, interpolation='nearest', cmap=cmap)#, aspect='auto')
plt.title(title, fontsize=20)
plt.colorbar()
if target_names is not None:
tick_marks = np.arange(len(target_names))
plt.xticks(tick_marks, target_names, rotation=45, fontsize=15)
plt.yticks(tick_marks, target_names, fontsize=15)
if normalize:
cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]
thresh = cm.max() / 1.5 if normalize else cm.max() / 2
for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):
if normalize:
plt.text(j, i, "{:0.4f}".format(cm[i, j]),
horizontalalignment="center",
color="white" if cm[i, j] > thresh else "black")
else:
plt.text(j, i, "{:,}".format(cm[i, j]),
horizontalalignment="center",
color="white" if cm[i, j] > thresh else "black")
plt.tight_layout()
plt.ylabel('True label', fontsize=15)
plt.xlabel('Predicted label\n misclass={:0.4f}'.format(misclass), fontsize=15)
plt.show()
X= data_clas.drop(['Retain','Profit0', 'District_first','AgeZero_0fascia','IncomeZero_0fascia'], axis=1)
y= data_clas['Retain']
x_train, x_test, y_train, y_test = train_test_split(X, y,
stratify=y,
test_size=0.3, random_state=2)
X_TRAIN= x_train
Y_TRAIN=y_train
Regressione Logistica¶
Per la regressione logistica, sono stati seguiti due diversi approcci:
Primo approccio:
- Si ricercano i pesi della regressione logistica tramite cross validation ( il valore ottimale dei pesi è stato scelto tramite f1-score, che bilancia precisione and recall).
- Si performa la regressione logistica su tutto il dataset train con i pesi ottenuti.
- Si calcola la matrice di confusione sul dataset test.
Secondo approccio
- Si performa una regressione logistica pesata, in cui i pesi sono fissati tramite il comando class_weight='balanced'
- Si sceglie il valore di soglia della probabilità ( oltre la quale assegnare 1 ) in due possibili modi, cioè viene scelto :
- il valore di probabilità che massimizza il G-mean ( una metrica che bilancia recall and specificità ) nel seguente modo: $$G-mean= \sqrt{recall \cdot\ specificità }.$$ Se i positivi sono classicati erroneamente (cioè sono predetti come negativi), il valore G-mean è basso anche se i negativi sono classificati correttamente.
- il valore di probabilità che massimizza F2-Measure ( una metrica che pone maggiore attenzione a minimizzare i falsi negativi ) , perchè è definita come $$ F2-Measure = \frac{5 \cdot Precision \cdot Recall}{4 \cdot Precision + Recall}$$.
- si calcola la matrice di confusione sul dataset test.
Regressione Logistica- I approccio¶
w = [ {0:1.0,1:1.0},
{0:1.0,1:0.1}, {0:10,1:0.1}, {0:100,1:0.1},{0: 1.0, 1: 5} , {0: 1.0, 1: 9},
{0:10,1:0.01}, {0:1.0,1:0.01}, {0:1.0,1:0.001}, {0:1.0,1:0.005},
{0:1.0,1:10},{0:1.0,1:50}, {0:1.0,1:99}, {0:1.0,1:100},
{0:1.0,1:200}, {0:1.0,1:1000}, {0:10,1:1000},{0:100,1:1000} ]
hyperparam_grid = {"class_weight": w}
lg3 = LogisticRegression(random_state=13)
grid = GridSearchCV(lg3,hyperparam_grid,scoring="f1", cv=100, n_jobs=-1, refit=True)
grid.fit(x_train,y_train)
print(f'Miglior valore di f1-score: {grid.best_score_} con parametri: {grid.best_params_}')
Performiamo la regressione logistica, con i pesi trovati dalla cross-validation e calcoliamo la matrice di confusione e alcune metriche sul dataset test.
lg3 = LogisticRegression(random_state=13, class_weight={0: 1.0, 1: 5})
#fit it
lg3.fit(x_train,y_train)
y_pred = lg3.predict(x_test)
plt.rcParams["axes.grid"] = False
lr_probs = lg3.predict_proba(x_test)
lr_probs = lr_probs[:, 1]
cm=confusion_matrix(y_test, y_pred)
target_names=['Clienti fedeli','Clienti persi']
plot_confusion_matrix(cm,
target_names,
title='Confusion matrix',
cmap=None,
normalize=False)
a=metrics(y_test, y_pred)
a
scores = pd.DataFrame(columns=['A'], index=['Regressione Logistica I','Regressione Logistica II','Regressione Logistica III','SVM','KNN'])
lr_probs = lg3.predict_proba(x_test)
lr_probs = lr_probs[:, 1]
forest_auc = roc_auc_score(y_test, lr_probs)
forest_fpr, forest_tpr, _ = roc_curve(y_test, lr_probs)
scores['A']['Regressione Logistica I'] = forest_fpr, forest_tpr, forest_auc
metriche = pd.DataFrame(columns=['Recall', 'Misclass'], index=['Regressione Logistica I','Regressione Logistica II','Regressione Logistica III','SVM','KNN'])
misclass=1- accuracy_score(y_test, y_pred)
recall = recall_score(y_test, y_pred)
metriche['Recall']['Regressione Logistica I']= recall
metriche['Misclass']['Regressione Logistica I']= misclass
from matplotlib import pyplot
from numpy import sqrt, argmax
model = LogisticRegression(solver='lbfgs', random_state=1, class_weight='balanced')
model.fit(x_train,y_train)
pyplot.rcParams["axes.grid"] = True
# predict probabilities
yhat = model.predict_proba(x_train)
# keep probabilities for the positive outcome only
yhat = yhat[:, 1]
# calculate roc curves
fpr, tpr, thresholds = roc_curve(y_train, yhat)
# calculate the g-mean for each threshold
gmeans = sqrt(tpr * (1-fpr))
# locate the index of the largest g-mean
ix = argmax(gmeans)
# plot the roc curve for the model
pyplot.plot([0,1], [0,1], linestyle='--', label='No-skill')
pyplot.plot(fpr, tpr, label='Logistic')
pyplot.scatter(fpr[ix], tpr[ix], marker='o', color='black', label='Best')
# axis labels
pyplot.title('ROC Curve', fontsize=20)
pyplot.xlabel('False Positive Rate', fontsize=20)
pyplot.ylabel('True Positive Rate',fontsize=20)
pyplot.legend(prop={'size': 20})
# show the plot
pyplot.show()
print('Il valore di probabilità che massimizza G-Mean è %f, con G-Mean=%.3f' % (thresholds[ix], gmeans[ix]))
model.fit(x_train,y_train)
y_pred = model.predict_proba(x_test)[:,1]
threshold= thresholds[ix]
y_pred = np.where(y_pred > threshold, 1, 0)
# Find prediction to the dataframe applying threshold
cm=confusion_matrix(y_test, y_pred)
target_names=['Clienti fedeli','Clienti persi']
plt.rcParams["axes.grid"] = False
plot_confusion_matrix(cm,
target_names,
title='Confusion matrix',
cmap=None,
normalize=False)
b=metrics(y_test, y_pred)
b
lr_probs = model.predict_proba(x_test)
lr_probs = lr_probs[:, 1]
forest_auc = roc_auc_score(y_test, lr_probs)
forest_fpr, forest_tpr, _ = roc_curve(y_test, lr_probs)
scores['A']['Regressione Logistica II'] = forest_fpr, forest_tpr, forest_auc
misclass=1- accuracy_score(y_test, y_pred)
recall = recall_score(y_test, y_pred)
metriche['Recall']['Regressione Logistica II']= recall
metriche['Misclass']['Regressione Logistica II']= misclass
(F2-Score)¶
Proviamo a cambiare soglia di probabilità .
model = LogisticRegression(solver='lbfgs', random_state=1, class_weight='balanced')
model.fit(x_train,y_train)
pyplot.rcParams["axes.grid"] = True
# predict probabilities
yhat = model.predict_proba(x_train)
# keep probabilities for the positive outcome only
yhat = yhat[:, 1]
# calculate pr-curve
precision, recall, thresholds = precision_recall_curve(y_train, yhat)
# convert to f score
fscore = (5 * precision * recall) / (4*precision + recall+ 0.000000000000000000001) #per non avere zero al denominatore
# locate the index of the largest f score
ix = argmax(fscore)
# plot the roc curve for the model
no_skill = len(y_train[y_train==1]) / len(y_train)
pyplot.plot([0,1], [no_skill,no_skill], linestyle='--', label='No-skill')
pyplot.plot(recall, precision, label='Regressione Logistica')
pyplot.scatter(recall[ix], precision[ix], marker='o', color='black', label='Best')
# axis labels
pyplot.title('Curva Precisione vs Recall', fontsize=20)
pyplot.xlabel('Recall', fontsize=20)
pyplot.ylabel('Precision', fontsize=20)
pyplot.legend(prop={'size': 20})
# show the plot
pyplot.show()
print('Il valore di probabilità che massimizza F2-score è %f, con F2-Score=%.3f' % (thresholds[ix], fscore[ix]))
model.fit(x_train,y_train)
y_pred = model.predict_proba(x_test)[:,1]
threshold= thresholds[ix]
y_pred = np.where(y_pred > threshold, 1, 0)
# Find prediction to the dataframe applying threshold
cm=confusion_matrix(y_test, y_pred)
target_names=['Clienti fedeli','Clienti persi']
plt.rcParams["axes.grid"] = False
plot_confusion_matrix(cm,
target_names,
title='Confusion matrix',
cmap=None,
normalize=False)
c=metrics(y_test, y_pred)
c
lr_probs = model.predict_proba(x_test)
lr_probs = lr_probs[:, 1]
forest_auc = roc_auc_score(y_test, lr_probs)
forest_fpr, forest_tpr, _ = roc_curve(y_test, lr_probs)
scores['A']['Regressione Logistica III'] = forest_fpr, forest_tpr, forest_auc
misclass=1- accuracy_score(y_test, y_pred)
recall = recall_score(y_test, y_pred)
metriche['Recall']['Regressione Logistica III']= recall
metriche['Misclass']['Regressione Logistica III']= misclass
SVM¶
Per i successivi algoritmi, normalizziamo il dataset, in quanto sono algoritmi che si basano su concetti di distanza.
scaler = StandardScaler()
a= np.asarray( x_train.columns.tolist())
b= np.asarray( x_test.columns.tolist())
X_TRAIN = scaler.fit_transform(X_TRAIN)
X_TRAIN=pd.DataFrame(X_TRAIN , columns=a )
x_test = scaler.fit_transform(x_test)
x_test=pd.DataFrame(x_test , columns=b)
Performiamo la cross-validation (K=5) per cercare i tuning parameters degli algoritmi.
In particolare, applichiamo una versione 'modificata' della cross-validation:
- in ogni iterazione, è esclusa una parte del dataset per la validazione
- le k-1 folds del training vengono bilanciate tramite undersampling ( RandomUnderSampler), mentre la fold Validation non viene modificata.
- Ripetiamo K volte
In questo modo, i dati in fase di validazione sono rappresentativi del dataset reale.
In sintesi, il modello è 'allenato' su un dataset bilanciato, mentre la validazione avviene su un dataset non bilanciato.
Dopo aver scelto il valore ottimale del parametro, si fitta il modello su un dataset bilanciato. Quindi si calcola la matrice di confusione sul dataset test non bilanciato.
Il valore ottimo del tuning parametro è scelto sulla base del valore di F2-score. In particolare:
- per SVM, cerchiamo il valore ottimo di C, che è un parametro di regolarizzazione. Come possibili valori di C, proviamo C=[0.1, 1, 10 , 100, 1000];
- per KNN, cerchiamo il valore ottimo di K, che è il numero di vicini che vengono considerati dall'algoritmo. Come possibili valori di K, proviamo K da 1 a 20.
from imblearn.under_sampling import RandomUnderSampler
spl=5
random_seed=1
kf = StratifiedKFold(n_splits=spl, random_state=random_seed)
C_range = [0.1, 1, 10, 100, 1000]
griglia= 5
cross_val_fbeta_score_lst= np.zeros(shape=(len(C_range),spl))
rus = RandomUnderSampler(random_state=0, replacement=True)
I boxplot che seguono mostrano le statistiche su F2-score ottenute durante la cross-validation per ogni valore di C. Come vediamo, il migliore valore di C è 10.
#svm
i=-1
for train_index_ls, validation_index_ls in kf.split(X_TRAIN, Y_TRAIN):
train, validation =X_TRAIN.iloc[train_index_ls], X_TRAIN.iloc[validation_index_ls]
target_train, target_val = Y_TRAIN.iloc[train_index_ls], Y_TRAIN.iloc[validation_index_ls]
rus = RandomUnderSampler(random_state=0, replacement=True)
X_train_res, y_train_res = rus.fit_resample(train, target_train)
i=i+1
for C in range( len( C_range)) :
model = SVC(C=C_range[C])
model.fit(X_train_res, y_train_res)
validation_preds= model.predict(validation)
cross_val_fbeta_score_lst[C][i]=fbeta_score(target_val, validation_preds, beta=2, pos_label=1)
df = pd.DataFrame(cross_val_fbeta_score_lst.transpose(),
columns= ['C=0.1','C=1', 'C=10', 'C=100', 'C=1000'])
#boxplot = df.boxplot( rot=45, fontsize=15)
#df.plot.box(return_type='axes', figsize=(12,8))
boxplot = df.boxplot(column=['C=0.1','C=1', 'C=10', 'C=100', 'C=1000'], fontsize=15, rot=45, figsize=(12,8))
# print ('Cross validated beta score: {}'.format(cross_val_fbeta_score_lst.mean(axis=1)))
vettore_medie=cross_val_fbeta_score_lst.mean(axis=1)
m=np.argmax(vettore_medie)
#cross_val_fbeta_score_lst
Quindi fittiamo il modello con C=1 su il dataset train, dopo averlo bilanciato.
x_train_rus, y_train_rus = rus.fit_resample(X_TRAIN, Y_TRAIN)
model = SVC(C=1, probability=True)
model.fit(x_train_rus, y_train_rus)
Riportiamo quindi le metriche e la matrice di confusione sul dataset test.
validation_preds= model.predict(x_test)
cm=confusion_matrix(y_test, validation_preds)
plt.rcParams["axes.grid"] = False
plot_confusion_matrix(cm,
target_names=['Clienti fedeli','Clienti persi'],
title='Confusion matrix',
cmap=None,
normalize=False)
d=metrics(y_test, validation_preds)
d
lr_probs = model.predict_proba(x_test)
lr_probs = lr_probs[:, 1]
forest_auc = roc_auc_score(y_test, lr_probs)
forest_fpr, forest_tpr, _ = roc_curve(y_test, lr_probs)
scores['A']['SVM'] = forest_fpr, forest_tpr, forest_auc
misclass=1- accuracy_score(y_test, validation_preds)
recall = recall_score(y_test, validation_preds)
metriche['Recall']['SVM']= recall
metriche['Misclass']['SVM']= misclass
KNN¶
Ricerchiamo il valore ottimo di K mediante cross-validation.
random_seed=1
spl=5
griglia=20
kf = StratifiedKFold(n_splits=spl, random_state=random_seed)
cross_val_fbeta_score_lst= np.zeros(shape=(griglia,spl))
#KNN
# Train the model using the training sets
i=-1
for train_index_ls, validation_index_ls in kf.split(X_TRAIN, Y_TRAIN):
# keeping validation set apart and oversampling in each iteration using smote
train, validation =X_TRAIN.iloc[train_index_ls], X_TRAIN.iloc[validation_index_ls]
target_train, target_val = Y_TRAIN.iloc[train_index_ls], Y_TRAIN.iloc[validation_index_ls]
rus = RandomUnderSampler(random_state=0, replacement=True)
X_train_res, y_train_res = rus.fit_resample(X_TRAIN, Y_TRAIN)
i=i+1
for k in range(1,griglia+1): #k=1 a k=griglia
model = KNeighborsClassifier(n_neighbors=k)
model.fit(X_train_res, y_train_res)
validation_preds= model.predict(validation)
cross_val_fbeta_score_lst[k-1][i]=fbeta_score(target_val, validation_preds, beta=2, pos_label=1)
#cross_val_fbeta_score_lst[k-1][i]=recall_score(target_val, validation_preds, pos_label=1)
#print ('Cross validated beta score: {}'.format(cross_val_fbeta_score_lst.mean(axis=1)))
vettore_medie=cross_val_fbeta_score_lst.mean(axis=1)
train_mean = np.mean(cross_val_fbeta_score_lst, axis=1)
train_std = np.std(cross_val_fbeta_score_lst, axis=1)
# Draw lines
plt.plot(range(1,griglia+1), train_mean, '--', color="#111111", label="Training score")
# Draw bands
plt.fill_between(range(1,griglia+1), train_mean - train_std, train_mean + train_std, color="#DDDDDD")
# Create plot
plt.title("cross validation Curve", fontsize=20)
plt.xlabel("K", fontsize=20), plt.ylabel("2-SCORE", fontsize=20), plt.legend(loc="best", prop={'size': 20})
plt.tight_layout()
m=np.argmax(vettore_medie)
k_ottimo=m+1 #nella posizione zero del vettore c'è k=1
plt.plot(k_ottimo,max(vettore_medie),'ro')
plt.xticks( range(1,griglia+1) ,range(1,griglia+1), fontsize=20)
plt.yticks(fontsize=20)
plt.show()
print('Il valore ottimale di K è 1.')
Alleniamo il modello sul dataset train bilanciato e riportiamo le metriche e la matrice di confusione sul dataset test.
x_train_rus, y_train_rus = rus.fit_resample(X_TRAIN, Y_TRAIN)
model = KNeighborsClassifier(n_neighbors=1)
model.fit(x_train_rus, y_train_rus)
# F0.5-measure puts more attention on minimizing false positives t
#positivo= 1
#voglio minimizzare chi è zero ma dice che è 1 , cioè falsi negativi
validation_preds= model.predict(x_test)
cm=confusion_matrix(y_test, validation_preds)
plt.rcParams["axes.grid"] = False
plot_confusion_matrix(cm,
target_names=['Clienti fedeli','Clienti persi'],
title='Confusion matrix',
cmap=None,
normalize=False)
e= metrics(y_test, validation_preds)
e
lr_probs = model.predict_proba(x_test)
lr_probs = lr_probs[:, 1]
forest_auc = roc_auc_score(y_test, lr_probs)
forest_fpr, forest_tpr, _ = roc_curve(y_test, lr_probs)
scores['A']['KNN'] = forest_fpr, forest_tpr, forest_auc
misclass=1- accuracy_score(y_test, validation_preds)
recall = recall_score(y_test, validation_preds)
metriche['Recall']['KNN']= recall
metriche['Misclass']['KNN']= misclass
Conclusioni¶
Abbiamo visto 5 diversi modelli di classificazione.
Il peggior modello è il KNN. I modelli di regressione logistica risultano essere, in generale, i migliori.
Per riassumere il lavoro di classificazione mostriamo in seguito il valore di misclass, ( 1- accuracy, cioè la percentuale di osservazioni mal classificate) , il valore di Recall ( che vale 1 se non ci sono falsi negativi) e la ROC Curve.
fig = go.Figure()
fig = fig.add_trace(go.Scatter(x=scores['A']['Regressione Logistica I'][0], y=scores['A']['Regressione Logistica I'][1], mode='lines',
line_shape='linear', name='Regressione Logistica I'))
fig = fig.add_trace(go.Scatter(x=scores['A']['Regressione Logistica II'][0], y=scores['A']['Regressione Logistica II'][1], mode='lines',
line_shape='linear', name='Regressione Logistica II'))
fig = fig.add_trace(go.Scatter(x=scores['A']['Regressione Logistica III'][0], y=scores['A']['Regressione Logistica III'][1], mode='lines',
line_shape='linear', name='Regressione Logistica III'))
fig = fig.add_trace(go.Scatter(x=scores['A']['SVM'][0], y=scores['A']['SVM'][1], mode='lines',
line_shape='linear', name='SVM'))
fig = fig.add_trace(go.Scatter(x=scores['A']['KNN'][0], y=scores['A']['KNN'][1], mode='lines',
line_shape='linear', name='KNN'))
fig.update_xaxes(range=[-0.05, 1.05], title_text='False Positive Rate')
fig.update_yaxes(title_text='True Positive Rate')
fig.update_layout(title='ROC Curve')
fig.show()
metriche